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Abstract 

A spline-based method for approximating thin shell dynamics is presented here. While the 
method is developed in the context of the Donnell-Mushtari thin shell equations, it can be 
easily extended to the Byrne- Fliigge-Lur’ye equations or other models for shells of revolution 
as warranted by applications. The primary requirements for the method include accuracy, 
flexibility and efficiency in smart material applications. To accomplish this, the method was 
designed to be flexible with regard to boundary conditions, material nonhomogeneities due 
to sensors and actuators, and inputs from smart material actuators such as piezoceramic 
patches. The accuracy of the method was also of primary concern, both to guarantee full 
resolution of structural dynamics and to facilitate the development of PDE-based controllers 
which ultimately require real-time implementation. Several numerical examples provide initial 
evidence demonstrating the efficacy of the method. 
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1 Introduction 


A spline-based method for numerically approximating thin shell dynamics is presented in 
this paper. This research is motivated by the growing number of applications involving shell 
models for which explicit solutions are typically unavailable. A specific example involves the 
development of control strategies for reducing vibrations, fatigue and structure-borne noise in 
aircraft fuselages. Typically, some form of cylindrical shell model is employed with the specific 
equations and boundary conditions dictated by the application or experimental setup. When 
numerically discretizing the models in such applications, the approximation method must be 
flexible with respect to boundary conditions and adaptable with regard to nonhomogeneities 
in materials and geometries. For example, if piezoceramic patches bonded to or embedded in 
the shell are used as actuators or sensors, the numerical method must be sufficiently flexible so 
as to permit extension to models incorporating these components. Also, the primary motion 
in many vibration and noise control applications is bending- dominated and the numerical 
method must accurately approximate such dynamics. Finally, from a control perspective, the 
approximation method must adequately preserve stability properties of the physical system 
as well as be sufficiently efficient so as to permit real-time implementation. 

Current techniques for approximating shell dynamics include modal expansions [17, 18, 
26], finite element discretizations [2, 3, 4, 11, 13, 14, 15, 20, 21, 22, 25, 28, 30] and finite 
difference approximations [31, 32]. From a theoretical perspective, modal expansions arise 
naturally when separating variables in linear shell models. In models for which analytic 
expressions for the eigenfunctions or modes can be obtained, this provides a straightforward 
and often quite efficient method of approximation. The difficulty, however, lies in the fact 
that analytic expressions describing mode shapes are known only for a restricted class of 
boundary conditions and for a limited set of models. For general boundary conditions and 
models, or systems involving coupling between a shell and an adjacent component (e.g., an 
acoustic field or a piezoceramic actuator), the modes must first be numerically approximated 
or experimentally determined before expansions can be formed. 

Experimental determination of natural frequencies and modes is typically accomplished 
by exciting the structure using an acoustic source, shaker, magnetic driver, et cetera. At each 
natural frequency, the shape of the corresponding mode is determined and characterized for 
use in modal expansions. As detailed in [32, 33], however, difficulties are encounted in this 
procedure in plate and shell-like structures due to the presence of multiple independent mode 
shapes which can occur at single experimental frequencies. In such cases, the experimentalist 
must excite the structure at various locations and complete an orthogonalization process to 
obtain a complete modal basis. The determination of which frequencies/modes to test for 
this behavior requires numerical analysis or extreme care when performing the experiments. 
Furthermore, as discussed in [32], internal or air damping will cause distortion of modes (e.g., 
modal lines that do not cross) which if not accounted for, will degrade modal expansions for 
approximating the structural dynamics. 

In applications, modes determined for boundary conditions yielding analytic expressions 
(e.g., simply-supported end conditions) are occasionally used to approximate solutions to mod- 
els with more general boundary conditions or models incorporating additional components. 
In some cases, the influence of various boundary conditions on modal characteristics has been 
investigated [17]; however, without any convergence analysis for such techniques, convergence 
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of the numerical approximates to the true model solution cannot be guaranteed. Finally, it 
is difficult to apply modal techniques in applications modeled by nonlinear shell equations 
or nonlinear actuator contributions (e.g., electrostrictive actuators or piezoceramic actuators 
embedded in constrained damping layers). 

Finite element approaches, on the other hand, are directly applicable for a variety of 
boundary and coupling conditions and in models including actuator contributions or nonlinear 
dynamics. Moreover, a large body of software exists for obtaining finite element solutions for 
various shell models (see [13, 20, 21, 22] for discussions of finite element methods and software 
for shells of revolution). 

A difficulty when developing and applying finite element methods for shell applications, 
however, concerns the manifestation of various locking phenomena. Shear locking , which 
has also been studied extensively in the context of Reissner-Mindlin plate models, is due to 
element incompatibility when enforcing the Kirchhoff-Love constraint of vanishing transverse 
shear strains as the shell thickness h tends to zero [1, 10]. In shell applications, an even more 
serious problem leading to the failure of various finite element methods is the phenomenon of 
membrane or shear-membrane locking. This phenomenon occurs when the total deformation 
energy is bending dominated, and is due to smoothness and asymptotic constraints in the shell 
model which are not appropriately represented by the approximation method (see [3, 4, 15, 25, 
28]). If these constraints are not satisfied by approximating elements, the numerical solution is 
often overly stiff in the sense that the model exhibits bending dynamics which the approximate 
solution cannot match. As detailed in [25], mesh sizes must be chosen significantly smaller 
than the shell thickness to ensure accurate approximations with high-order finite elements in 
such bending dominated regimes. These examples also illustrate that even with such mesh 
size restrictions, low-order finite element methods often fail in such regimes. The use of finite 
elements which exhibit locking is detrimental in simulations and will seriously degrade the 
performance of any controllers designed around such methods. 

A second issue which must be considered when developing finite element methods for shells 
concerns the piecewise constant material parameters and inexact boundary conditions which 
often arise in smart material applications. For example, the use of piezoceramic actuators leads 
to piecewise discontinuities in the density, stiffness, Poisson ratio and damping parameters. 
Consequently, finite element meshes must be aligned with the actuator boundaries to maintain 
accuracy (the reader is referred to [34] for a discussion of finite elements designed specifically 
for piezoelectric applications). In experiments involving shells, slight energy loss must often be 
modeled into the boundary conditions to account for inexact boundary clamps. Consequently, 
elements must be provided with extra degrees of freedom to accommodate these boundary 
conditions. Standard commercial codes not providing these capabilities will lead to potentially 
inaccurate results when applied to shells with smart material actuators and/or sensors. 

The use of finite difference discretizations for approximating thin shell dynamics is less 
common due to inherent difficulties with the high-order equations and boundary conditions 
arising in the models. The reader is referred to [31] for further discussion of finite difference 
methods for shell applications. 

In this paper, we consider a Galerkin method for discretizing thin shell models with linear 
or cubic splines chosen as basis elements for approximating the longitudinal and circumferen- 
tial motion and cubic splines used to approximate the transverse component. These choices 
are motivated by smoothness and accuracy criteria, adaptability for a variety of boundary 


2 



conditions, and suitability when including the contributions of actuators such as piezoceramic 
patches bonded to the shell [6]. For the discussion here, the shells are assumed to have lengths 
that are relatively short in relation to the radii, and the Donnell-Mushtari shell equations are 
employed. This facilitates the discussion while providing a framework which is easily extended 
to more accurate long shell models through the inclusion of Byrne- Fliigge-Lur ye components. 
In this study, the formulation of the method and numerical examples demonstrating the ac- 
curacy, efficiency and flexibility of the method are presented. Emphasis in these examples 
is placed on demonstrating that when h x denotes the axial mesh size, the expected 0(h “.) 
and G{h*) accuracy of the method is maintained when approximating with linear and cubic 
splines, respectively. Convergence analysis and further analysis of the method with regards 
to membrane and shear locking will appear in a future companion paper. 

In Section 2, the strong and weak forms of the modeling equations for a thin shell with 
surface-mounted piezoceramic patches are summarized. Modal analysis for the special case 
of an undamped shell with constant coefficients and simply-supported boundary conditions is 
discussed in Section 3. This provides a framework for testing the convergence of the approxi- 
mate mass and stiffness components in the system. The approximation method and resulting 
finite dimensional matrix system are detailed in Section 4. Finally, examples illustrating the 
method are presented in Section 5. These include modal approximations for various bound- 
ary conditions as well as results demonstrating convergence rates when approximating shell 
dynamics generated by various external inputs. To illustrate the shell dynamics generated by 
the piezoceramic patches, a modeled voltage spike to the patches is used as input in the final 
example. The resulting frequencies are then compared with those obtained in the eigenvalue 
(modal) analysis to demonstrate the consistency of the method. 


2 Donnell-Mushtari Equations 



X 

Figure 1 . Thin cylindrical shell with surface mounted piezoceramic patches. 
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We consider here a shell of length £, thickness k and radius R having mass density p, 
Young’s modulus E, Poisson ratio v and Kelvin- Voigt damping coefficient c D . As depicted 
in Figure 1, the axial direction is taken along the x-axis. The displacements of the middle 
surface in the longitudinal, circumferential and transverse directions are denoted by u, v and 
w, respectively. 

Bonded to the shell are s pairs of piezoceramic patches which can be employed as sensors 
and/or actuators in structural and structural acoustic applications [5, 19]. To simplify the 
exposition, the patches are all assumed to have thickness h pe , Young’s modulus E pe , Poisson 
ratio v pe , and Kelvin- Voigt damping coefficient c Dpe . Furthermore, it is assumed that the 
glue bonding layer provides negligible contribution to the structural dynamics. The reader is 
referred to to [6, 9] for details concerning the incorporation of differing patch characteristics 
and bonding layers in the ensuing models). 

Throughout this discussion, it will be assumed that external inputs to the shell will be in 
the form of transverse, longitudinal and circumferential surface forces as well as line moments 
and forces generated by the patches. The surface forces can be used to model a variety of 
phenomena including coupling interactions with adjacent fields (e.g., acoustic fields - see [5]) 
and input from certain actuators. The patch moments and forces arise when the elements are 
used as actuators. While more general inputs can be considered, the above-mentioned cases 
demonstrate the flexibility of the numerical method for typical smart materials applications. 


2.1 Strong Form of Equations 


As detailed in [6, 26], moment and force balancing yields the Donnell-Mushtari equations 
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as a model for the thin shell dynamics. Here M x , Mg, M x $, and Mg x are internal moments, 
while N x , Ng, N x g and Ng x denote internal force resultants. External surface forces are denoted 
by q x , qe, Qn whereas the external resultants (line moments and forces) generated by the i th 
patch pair are designated by (M x ) pei ,(M 9 ) pei ,(N x ) pei and (7V*) pei . The indicator function 

S pei (x,6) = Si, 2 (x)S l!2 (0) , 

where 
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indicates the sense of the forces generated by the i th patch pair having edges at , x 2 »» 0 \ i, 02 i ■ 

The symmetry of the function arises from the property that for homogeneous patches hav- 
ing uniform thickness, equal but opposite strains are generated about the point = 

((xi, + x 2 ,)/2, ($ 1 ,- + 6 - 2 i)/ 2 ) in the two coordinate directions. Similarly, the characteristic 
function 

1 , xu < x < x 2 , , Ou < 9 < 02 i 
0 , otherwise 

will be used in ensuing discussion to isolate internal and external contributions due to the i th 
patch pair. 

Under the assumption that stress is proportional to a linear combination of strain and 
strain rate, the internal moments and force resultants in regions of the shell not covered by 
patches are given by 
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(see [6] or [26] for the undamped case). The substitution of these equations in (2.1) yields 
the damped form of the Donnell-Mushtari equations for a uniform, homogeneous shell that is 
devoid of patches. 

The bonding of piezoceramic patches to the shell produces contributions due to both 
internal (material) and external moments and forces. The internal contributions are due to 
the geometrical and material changes afforded by the patches. The external contributions 
result from the piezoelectric effect in the patches which is manifested as generated strains in 
response to applied voltages. 

As detailed in [6, 9], the internal force resultant N x is given by 
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for systems in which the patches have identical properties (compare with (2.2)). Analogous 
expressions are given in [6, 9] for the remaining resultants and the resultants when the patches 
have differing material characteristics. 

To characterize the external contributions, it is typical to start with the assumption that 
the strains generated by a patch are proportional to the applied voltage [6]. Since differing 
voltages can be applied to the outer and inner patches in the pair, we will differentiate between 
the two with VJi(t) and V, 2 (t) used to denote the voltages to the outer and inner patches in 
the i th pair, respectively. The proportionality constant relating generated strain to the input 
voltage is designated by d-$\. As detailed in [6], the total external moments and forces generated 
by the patches are 


[M x ) pei = [(M x ) peii + [M x ) pea ] Xpe ,(x,0) 

( M *)pe, = [( M s) petl + ( M e)pe i2 \ Xpe i{x,9) 

Wpe, = [(^pe* + WpeJ Xpe, (*, 0)S pei {x, 6) 
Wpe, = [C^Ope.-, + WpeJ Xpe i {x,0)S pei (x, 9) 

where 
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When substituted into (2.1), the expressions (2.4) provide the input from the patches when 
voltages are applied. 
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2.2 Boundary Conditions 

Appropriate boundary conditions are dictated by the experimental setup or application under 
consideration. In many setups, such as the experimental shell apparatus at NASA Langley 
Research Center, the shell is supported by heavy endcaps. In such cases, the fixed-edge 
conditions 

dw 

u = v = w = — - = 0 , x = 0,£ ( 2 - 6 ) 

Ox 

may adequately model the end behavior of the shell. 

If slight boundary movement or rotations are suspected, “almost fixed” boundary condi- 
tions of the type discussed in [7, 8] can be employed. Such boundary conditions can be used 
to model the slight energy which result due to imperfect clamping of the structure. 

A third type of commonly considered boundary conditions are the simply-supported or 
shear diaphragm edge conditions 

v = w = M x = N x = 0 , x = 0,1 . (2-7) 

These boundary conditions are theoretically attractive since they provide one of the few sit- 
uations in which explicit modal expansions can be calculated for the Donnell-Mushtari shell 
equations. They are of limited use when modeling experimental shells, however, since they 
are appropriate only for endcaps which prevent movement in the v and w directions but 
are sufficiently flexible in the x-direction so that negligible moments M x and forces N x are 
generated. 

The Galerkin method of this work is equally efficient to implement for models incorporating 
the fixed-edge boundary conditions (2.6), “almost fixed” boundary conditions, or the simple- 
supported edge conditions (2.7). The method is demonstrated for the fixed-edge and simply- 
supported conditions while modifications to adapt the method to “almost fixed” can be found 
in [7, 8]. Later discussion will also illustrate the manner through which the method can be 
adapted to alternative boundary conditions which many arise in physical applications. 

2.3 Weak Form of Equations 

As noted in (2.1), the strong form of the modeling equations involves first and second deriva- 
tives of both the moment and force resultants. For structures with surface- mounted or embed- 
ded actuators or sensors (e.g., piezoceramic patches), this differentiation leads to difficulties 
due to discontinuities in the resultants. The internal resultants contain piecewise discontinu- 
ities due to material changes introduced by the patches (see for example, (2.3)). The external 
patch contributions are also discontinuous since they are defined only over regions covered by 
patches (see (2.4)). As a result, formal analysis and approximation using the strong form of 
the modeling equations leads to derivatives of the Dirac delta “function.” 

To alleviate these difficulties, we consider a weak form of the modeling equations. As 
detailed in [6], such equations can be derived from energy considerations and are equivalent to 
the strong form of the modeling equations under suitable smoothness conditions. This yields 
a form of the model which facilitates well-posedness analysis [5] and eliminates the difficulties 
associated with the discontinuous resultants. Moreover, smoothness requirements on basis 
functions are reduced which proves advantageous when constructing system matrices. 
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The weak form of (2.1), as derived from energy considerations in [6], is given by 

I L { Rph W’ h + RN ’l^ + N ” W - ~ = 0 

L l + "•w + RN “ji ~ K > ,r > 2 - gw**.§j } dxde = 0 
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for all ff = [t}i,t) 2 ,t} 3 ] £ V where V denotes the space of test functions (as detailed in [6], the 
indicator functions appear only in the definition of the external force resultants in the weak 
form of the equations). Specific choices for V depend upon the boundary conditions. For the 
fixed-edge boundary conditions (2.6), the space of test functions is taken to be 


where 


v = Him) x Him) x Him) 

= {v e H\si ) : ,(o) = nW = 0} 

Him = {■? 6 H 2 ((l) : 1,(0) = -,,(0) = ,(£) = ,,(*) = o} . 

In the case of the simply-supported edge conditions (2.7), essential boundary conditions are 
imposed only on the circumferential and transverse functions and V is taken to be 


V = H\Q) x ffj(ft) x Hl(Sl) 


with 

Him = b e H 2 m ■ i(o) = vw = 0} . 

In general, V is simply taken as the subset of the traditional Sobolev spaces satisfying essential 
boundary conditions (the reader is referred to [7] for modifications to account for “almost 
fixed” boundary conditions). 

A comparison between (2.8) and (2.1) illustrates that in the weak form, derivatives are 
transferred from the discontinuous resultants onto suitably smooth test functions. This allevi- 
ates the difficulties associated with the discontinuities and reduces smoothness requirements 
on approximate solutions. 
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3 Modal Solution — Simply-Supported Boundary Con- 
ditions, Constant Coefficients, No Damping 


As noted in the introduction, closed form modal expressions for shell models can be obtained 
only in a limited number of cases. One case in which analytic frequencies and modes can 
be calculated is the undamped (cp = 0) shell model with simply-supported edge conditions 
(2.7) and constant parameters p, v and E. Separation of variables in this setting is classical 
and can be found in numerous texts (e.g., [27]); we summarize the arguments here to facili- 
tate numerical examples in later sections. While such conditions are not attained in typical 
shell applications, the consideration of the undamped, constant coefficient shell model with 
simply-supported edge conditions provides an excellent means of testing discretization tech- 
niques since approximate solutions can be compared with analytic values. The discretization 
techniques can then be used to approximate natural frequencies, modes and shell dynamics in 
models which incorporate piecewise constant parameters (including damping) and physically 
realistic boundary conditions. 

Throughout this section, it is assumed that cp = 0 and p, v and E are constant. Since 
our interest here is restricted to the calculation of natural frequencies and modes, we will also 
consider the shell model to be unforced (no external patch contributions and q x = qe = q n = 0). 
In this case, the strong form of the Donnell- Mushtari thin shell equations (2.1) can be written 
in the operator format 

ph— = Lu (3.1) 

where u = [u, v, w] T . The operator L here is given by 


8 2 (1 - u) d 2 

dx 2 + 2 R 2 dd 2 

Eh (1 + v) d 2 

1 - v 2 2 R dxdd 

v d 

-Rdi 


(1 + y) d 2 

2 R dxd6 
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R 2 d6 R 2 


v d 

Rdi 

\_d_ 

R? de 

I'd 4 2d 4 Id 4 
[dx 4 + R 2 dx 2 d6 2 + R* 


where k = h 2 / 12. Boundary conditions for the shell are denoted by 

Bu = 0 . 

Due to linearity, the displacements are expressed as 

u(t,x,6) = T(t)U(x,0 ) 
v(t, x,0) = T(t)V(x,0) 
w{t,x,d) = T(t)W{x,9 ) , 


and spatial and temporal components are separated to yield the eigenvalue problem 

LU + phu 2 U = 0 
BL 1 = 0 


(3.2) 
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along_with the temporal equation T" + u> 2 T = 0. Here u is the circular frequency of vibration 
and U = [U, V, WY contains the vibration modes in the axial, circumferential and transverse 
directions. We point out that separating variables to arrive at the eigenvalue problem (3.2) is 
equivalent to the assumption of a harmonic response in all components and one obtains the 
same eigenvalue problem in both cases. 

The structure of the eigenfunctions U,V and W is dependent upon the boundary conditions 
with closed form expressions attainable in only a few cases. To illustrate, we consider the 
eigenvalue problem (3.2) with the simply-supported boundary conditions (2.7). One form of 
the natural vibration modes for this case is 

U(x,6 ) = A\ cos ^ cos(m0) 

V ( x , 9) = Bi sin sin (m6) (3.3) 

W(x , 6) = C\ sin ^ ^ cos(m0) . 

The choice of the cosine representation for the x component of the axial vibration and sines 
for the circumferential and radial vibrations guarantees that the modes satisfy the boundary 
conditions. The relationship between the circumferential and axial/radial vibrations dictates 
that the former must be 90° out-of-phase from the latter two in 6 (see the moment and force 
expressions (2.2) or operator definition (3.1)). 

To determine the relationship between the frequencies u> and the wave numbers m and n. 
the expressions (3.3) are substituted into (3.2) to yield the system 

fl 2 — ^r(l + 1 ') v\ A\ 0 

TK1 + *) &-H 2 -m Bi = 0 (3.4) 

v\ — m Cl 2 — H 3 Ci 0 

where 

H i = A 2 + ^(1-„) 

H 2 = —(1 - is) + m 2 
H 3 = 1 + k (A 2 + m 2 ) 2 

and 

n 2 = 4(i - u 2 ) u 2 r 2 
t/ 

x mrR 
= ~L~ 
r k__ h 2 
~ R 2 ~ 12 R 2 ' 

A nontrivial solution is determined by setting the determinant to zero. This yields the 
cubic equation 

n 6 - K 2 n 4 + Kitf - K 0 = 0 ( 3 . 5 ) 
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in Q 2 . The coefficients here are given by 

K 2 = 1 + ^(3 - v) (m 2 + A 2 ) 4 - k ( m 2 + A 2 ) 2 

K x = ^(1 - v ) (3 + 2z/)A 2 + m 2 + (m 2 + A 2 ) 2 + j — ^-k ( m 2 + A 2 ) 3 j 

/To = — (1 — v) (1 — i/ 2 )A 4 + k (m 2 + A 2 ) j . 

The solutions to this cubic equation are then given by 


n 2 mnj = --y/K 2 -3Kx cos 


a + 2 (j - 1 )tt \ K 2 
3 / + 3 


where 


a = cos 


_j -27K 0 -2KI + 9KxK 2 


\ 2y/(Kl-ZKx) 3 ) 

Finally, the natural frequencies for the shell, with units of hertz, are given by 


f “• mnj 

Jmnj ~ ~ n 


, j = 1,2,3. 


2-kR V p{\ - v 2 ) ’ " 

We point out that due to the cubic nature of the characteristic equation, three natural fre- 
quencies and mode shapes are obtained for each set of wave numbers m and n. This leads to 
a significant interlacing of frequencies. 

A second set of modes which satisfy the boundary conditions and 0 compatibility criteria 


( TtTTX \ 

■ - ) sin(m$) 

( TL7TX \ 

— — J cos (m0) 


W(x,0) = C 2 sin — J sin (m0) . 

The system which arises when using this form of the modes is 

- n 2 -Hx -^(i + v) v\ i r a 2 
-^(l + r) W-H 2 m B 2 

i/A m Cl 2 - H 3 C 2 


which differs in sign in the (1,2), (2, 1), (2, 3) and (3, 2) elements from (3.4) which was obtained 
using the first set of modal expansions. The characteristic equations, however, are the same 
in both cases and are given by (3.5). Hence both sets of expressions yield the same natural 
frequencies for the shell. This should be expected since the second set simply represents a 
phase shift in 9 of the first. For this reason, most authors consider only the first set of vibration 
modes. 

In some instances, however, the second set (3.7) contributes linearly independent modes 
and hence must be retained in order to obtain a complete basis for approximation. This is 
illustrated in the examples of Section 5.1. Further details regarding the modal analysis can 
be found in [16]. 
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4 Approximation Method 

To approximate the solutions u,v,w to the system (2.8) or eigenvalue problem (3.2), basis 
functions B Uk (0, x), B Vk (0, x) and B Wk (0,x) in V are chosen and used to form expansions 

u N (t,e,x) = J2Mt)Bu k (0,x) 

k= 1 
-Vv 

v N (t,e,x) = Y,v k (t)B Vk (0,x) (4.1) 

k= 1 

w N (t,e,x)=Y i w k (t)B Wk (6,x). 

k = 1 

The time-dependent generalized Fourier coefficients Ufc(t), ujt(f), wjt(t) are then determined by 
orthogonalizing the residual with respect to linearly independent test functions from the span 
of the bases. 

To exploit the tensor nature of the cylindrical shell domain fi and periodicity in 0, the 
bases for the three displacements are constructed with Fourier components in 0 and spline 
components in x. When expressed with real components, this yields 

Nu 

U N (t,d,x) = u 0 n(t)B Un (x) 

n — 1 

Mu Nu Mu Nu ( 4 ' 2 ' 

+ ^2 u mn(t) cos(m6)B Un (x) + ^ ^2 u mn (t)sm(m0)B Un (x) 

771 = 1 71=1 771 = 1 71=1 

with similar expressions for v N (t, 0, x) and w N (t, 0, x). In the expansions, M U ,M V , M w are the 
Fourier limits and N u , N v , N w denote the number of splines B Un (x), B Vn (x), B Wn (x) used when 
approximating the longitudinal, circumferential and transverse displacements, respectivelv. 
The total number of basis functions for each component is Af u = N u - (2M U + 1),AT V = N„- 
(2M v + l) and J\f w = N w ■ (2 M w + 1). 

4.1 Axial Bcisis Functions 

The choice of splines in the axial variable is motivated by the following criteria: 

(i) efficiency; 

(ii) flexibility with regard to internal and external patch contributions; 

(iii) adaptability with respect to various boundary conditions; 

(iv) accuracy. 

For the longitudinal and circumferential displacements, both linear and cubic spline bases are 
considered. Due to differentiability requirements, only cubic splines are used when approxi- 
mating the transverse displacement. 

In all cases, a uniform partition along the x-axis is considered with gridpoints x n = 

nh x , h x = ij N ,n = 1, • • • , N. For n = — 1,0, 1, • • • , TV + 1, standard cubic splines are de- 
fined by 
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M x ) 



(x - x n _ 2 ) 3 , 

h 3 x + 3 hl(x - x n -i) + 3 h x (x - x„_i) 2 - 3(x - x n _ 2 ) 3 , 
h x “{“ 3 ^(Xn-)-] x) ~f~ 3 ^x(*^n+l X ) 3 (x 71 _j_i x) , 

(x n+2 - x) 3 , 


l 0 


(see [29]). The standard linear splines are defined by 


M x ) 



(x X n _i) 
(Xn-fl x) 
0 


•y x G [x n _i,x n ] 
, x G [x n , x n ^_i] 
, otherwise 




£ [^n — 2? l] 

^ £ [^n— 1 ^ -^n] 

G [^n? *£n+l] 

x G [x n+1 ,x n+2 ] 
otherwise 

(4.3) 

(4.4) 


for n = 0, * * • , N. 

In both cases, the standard splines must be modified to satisfy the essential boundary 
conditions to ensure approximate solutions in V . These modifications are summarized as 
follows. 


(i) Fixed displacements at x = 0,£: 


To construct basis functions satisfying fixed displacements (but unspecified slopes) at 
x = 0, the modified cubic splines are taken to be 

( &o(x) - 46_i(x) 

6i(x) - 6_i(x) 

^n(x) = ^ ^n(x) 

6^_i(x) — 6jv+i 
tjv(x) — 46/v + i(x) 


n = 0 
n = 1 

7i = 2, • • • , N — 2 

7 Z = iV — 1 

n = N 


(4.5) 


for a total of N + 1 functions. To satisfy the same condition, the first and last linear splines 
are omitted from the set to yield 

Cn(x) = c n (x) , n = l,---,iV- 1 . ( 4 . 6 ) 

It should be noted that with these definitions, 6„(0) = b n {£) = 0 and c„(0) = c n (£) = 0. 


(ii) Fixed displacements and slope at x = 0,£: 

Only the cubic splines are required to satisfy a condition of fixed displacement and slope 
at the endpoints since this is a condition imposed only on the transverse displacements. In 
this case, the modified splines are taken to be 

{ b 0 (x) — 26_i(x) — 2£>i(x) , 7i = l 

b n (x) , n = 2, - • ■ ,N — 2 , (4.7) 

&at(x) — 26;v_ 1 (x) — 26/v+i(x) , n = N - 1 

for a total of N — 1 basis functions. Note that these functions satisfy 

MO) = MO) = b n (£) = b' n (£) = 0 . 
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4.2 Bases for Specific Boundary Conditions 

Appropriate bases for various boundary conditions are then constructed by considering mod- 
ified splines satisfying any essential boundary conditions. For example, the bases used for 
simply-supported shells (2.7) must satisfy the conditions £ Vn (0) = B Vn (£) = 0, = 

= 0 (^e moment and shear conditions are natural and hence do not need to be explic- 
itly enforced). Choices for shell models with fixed-edge conditions are summarized in Table 1 
while corresponding choices for simply-supported edge conditions are given in Table 2. The 
last column in each table summarizes the total number of axial functions in each expansion. 

Once the axial bases are chosen, they are combined with the Fourier components in 0 
to form the total bases. We reiterate that the basis limits in each case are given by M u = 
N u ■ (2 M u + 1) = N v ■ (2 M v + 1) and = N w ■ (2 M w + 1). 

We point out that the linear and cubic spline bases described here are but two choices from 
among many that can be made for the axial components. For the applications considered here, 
cubic splines provided a good balance between accuracy, efficiency and adaptability with re- 
gards to patches and boundary conditions. If higher accuracy is desired, however (with slightly 
more expense when constructing system matrices), quintic splines can be employed. Similarly, 
spectral expansions employing Legendre, Chebyshev or sine functions can be employed once 
modifications have been made for boundary conditions. 


Shell Component 

Axial Basis Functions 

Component Definition 

Axial Limit 

longitudinal - linear 

= ^n(^) 

Cn(x) defined in (4.6) 

II 

£ 

1 

cubic 

B Un (x) = b n (x ) 

b n (x ) defined in (4.5) 

N u — N u + 1 

circumferential - linear 

B v n (x) = Cn(x) 

Cn(x ) defined in (4.6) 

r-H 

1 

fe? 

II 

cubic 

II 

TT 

of 

b n (x) defined in (4.5) 

Ny Ny + 1 

transverse - cubic 

3 

IT 

II 

o-i 

3 

IT 

6 n (x) defined in (4.7) 

1 

fe- 

ll 

<fe? 


Table 1. Axial basis definitions for models with fixed-edge boundary conditions (2.6). 


Shell Component 

Axial Basis Functions 

Component Definition 

Axial Limit 

longitudinal - linear 


Cn(x) defined in (4.4) 

N u = + 1 

cubic 

B u „(x) = bn(x) 

b n (x ) defined in (4.3) 

N u ~ N u + 3 

circumferential - linear 

Bv n (x) = Cn(x) 

Cn(x) defined in (4.6) 

1 

fe; 

ii 

<fef 

cubic 

B Vn (x) = b n (x) 

b n (x ) defined in (4.5) 

N V = N V + 1 

transverse - cubic 

= &n(^) 

S n (ar) defined in (4.5) 

T“ 1 
+ 

fe- 

ll 

'fe; 


Table 2. Axial basis definitions for models with simply-supported boundary conditions (2.7). 
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4.3 Components in the Finite Dimensional System 

With bases thus defined, approximating subspaces are taken to be = span{f? Ufc }£L , = 

spanjB^j^j and H ^ = spanjB^}^. It should be noted that H N = H ^ x x H £ C V. 
The approximating system is then determined by restricting the weak form (2.8) to H N with 
basis functions used as test functions. This is equivalent to orthogonalizing the residual with 
respect to elements from H N . 


4.4 Matrix System 


To form the resulting matrix system, the generalized Fourier coefficients for the three expan- 
sions (4.1) are first consolidated in the vectors 



Ui(t) 


MO 


MO 

U^{t) = 

j 

II 

to 

* 

: 

CH- 

II 









The full set of coefficients is then represented as 1 9^(t) — V jV "(<), W jVw (i)] r , where 

Af = Af u + Af v + Af w . 

The mass, stiffness and damping matrices as well as the forcing vector for the full system 
have the form 



’ U M 



M m = 


V m 





W m _ 


and 



Un + U12 

Vn 

+ v 12 

W 11 

Ki = 

U -21 + U22 

Vn 

+ V22 

W 21 


u 3l 

Vn 

EL 1 w 3k _ 



’ U u + U 12 

V u 

+ v 12 

Wn 

Kf = 

CD 

U2I + U22 

Vn 

+ v 22 

W21 


U31 

Vn 

ELi W» 


F u 


F*(t) = 


F v 

F w 


(4.9) 


(4.10) 


The various submatrices contain individual components which arise when the weak form (2.8) 
is restricted to H N . For example, the approximation of the mass and stiffness components in 
the longitudinal equation of (2.8) yields 
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(i) 


n t J 

ph + 2 y ] PpehpeXpei 


n l 

. 


EhR J^2E pe h pe R 

1 _ 1,2+22 l _ l/ 2 Xpe, 

t=l P e 


RB Uk B Uj dxd9 

dB Uk dB U] 


/•••\ r T r l r2ir rt J" Jjjfll, 2 E pe H pc U pe 

(m) 


5a: 5a; 

5ff^ 5g tij 
50 5a; 


dxdO 


dxdO 


(iv) WJa-jf + 


dB u 

B Wk ~^dxde 


, . n , . f* Eh ^ 1 

(v) pr»U = J J 5nT^ + Si 


Eh A E pe h pe 
i = i x + ^ 


2(1 + t') 


Xpe, 


5x 50 




Epehp e dB Uk dB Uj j 

[2/2(1 + v) + E R(1 + Vpe ) x *«\ 36 d9 dxde 


f2ir rl f * M "I 

(vii) [F u ]j = J o j o ^Rq x B Uj + Rj2{^x) pe ,-Q^->dxd6 


with similar expressions for the submatrices f/ n , V u , W u , V 12 and U l2 which contain the in- 
ternal damping contributions. The submatrices for the approximated circumferential and 
transverse equation are constructed analogously with details given in [16]. It should be noted 
that due to the presence of the characteristic functions in the force and moment resultants, 
the domains of integration for the patch contributions are restricted to those regions covered 
by patches. This is consistent with the fact that material and geometric changes, as well as 
patch inputs, occur only in those regions. 

In first-order form, the system has the form 


Kg 

0 

' d*{t) ' 


0 

0 


. &*{*) . 


-Ki 

' Kg 

0 

' t^(0) ' 


'aT 

0 

M* 

_ tK(0) _ 


. 92 . 


-K 


x 9*(t) 
x 9*(t) 


+ 


0 

L F"(t) 


(4.11) 


Multiplication of the inverted mass matrix yields a Cauchy equation of the form 

y N (t) = A N y N (t)+g N (t) 

y N ( o) = yff, 


(4.12) 
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where y N £ R 2J ^ ^(i)] • This yields a form which is suitable for simulations, 

parameter estimation and control applications. Note that the system can also be adapted to 
alternative boundary conditions through modifications of the first and last basis functions. 
Flexibility in this regard is also a hallmark of the method. 

Finally, we point out that approximation of the eigenvalue problem (3.2) using the Fourier- 
Galerkin method yields the generalized matrix eigenvalue problem 

Kgd* = u 2 M* r ti* r (4.13) 

where Kj[ and M are defined in (4.9). As illustrated in the examples of Section 5.1, (4.13) 
can be used to approximate the frequencies and modes for undamped shells with various 
internal characteristics (due to actuators and sensors) and boundary conditions. 


5 Examples 

To illustrate various facets of the Fourier- Galerkin method, two types of examples are pre- 
sented. In the first set of examples, approximation of frequencies and modes is considered 
for simply-supported and fixed-edge shells. This is accomplished by forming the mass and 
stiffness matrices and solving the generalized matrix eigenvalue problem (4.13). This serves 
two purposes. It illustrates the accuracy of the method in approximating the mass and stiff- 
ness components in the system and demonstrates the method as an effective technique for 
obtaining modal information when analytic or experimental values are unavailable. 

The second set of examples illustrates the approximation of static and dynamic shell 
responses to various forces. The initial steady state example illustrates that the expected 
G{h? x ) and 0(h x ) convergence rates are obtained when approximating a static shell response 
(this example involves the solution of a matrix system with coefficient matrix K g and vector 
F^). Analogous dynamic examples demonstrate that the same accuracy is obtained when the 
method is used to obtain the time-dependent system (4.12) which is then marched in time. 
The latter example also incorporates the internal Kelvin- Voigt damping. The dynamics of the 
system in response to patch inputs are considered in the final example. Natural frequencies are 
compared to those obtained by solving the generalized matrix eigenvalue problem to illustrate 
consistency among techniques. 

In these examples, the following shell characteristics were used: thickness h = .01 in , 
length i = 12 in, radius R = 3.0 in, mass density p = .283 lb/in 3 , Young’s modulus E = 3.0 x 
10 7 lb/in 2 , Poisson ratio v — 0.3 and Kelvin- Voigt damping coefficient cp = 15.09936 lb-in-s 2 . 
These values were chosen so as to permit comparison with Donnell-Mushtari shell results in 
[24, pp. 307-310]. We emphasize that the method is not restricted to such dimensions or 
ratios, however, and analogous convergence results have been obtained with a wide range of 
parameters. 

5.1 Modal Examples: Simply-Supported Shell 

In Section 3, separation of variables was used to obtain analytic expressions for modes and 
natural frequencies of shell models with simply-supported boundary conditions. Here we com- 
pare the approximate solutions obtained by solving the matrix eigenvalue problem (4.13) with 
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the analytic values. This provides a means of testing the accuracy of the method before using 
it to approximate in settings in which analytic solutions are unavailable. To illustrate, we 
consider the cases of purely axisymmetric modes (m = 0,n > 1), purely extensional modes 
(m > l,n = 0) and general modes (m > l,n > 1). 

5.1.1 m = 0,n > 1 Axisymmetric Modes 

Analytic frequencies for the axisymmetric case can be calculated via (3.6) with m = 0 while 
corresponding modes can be determined from (3.3) or (3.7) (see [16] for details). Corre- 
sponding approximate frequencies, obtained by solving (4.13) with Fourier limit M - 0, are 
compared with analytic values in Table 3. As noted in Section 3, three frequencies are ob- 
tained for each configuration of 1/2 wavelength (fixed values of n), and this is reflected in the 
analytic and approximate results in Table 3. Two of the frequencies correspond to coupled 
axial/radial modes while the remaining one corresponds to an uncoupled torsional (circum- 
ferential) mode. The frequencies of the torsional modes are boxed to facilitate comparison 
between the corresponding functional modal approximations. 

The cubic spline approximates in Table 3 were obtained using N - 8 basis functions 
whereas N = 32 linear splines were used to obtain the corresponding results in columns 8-10. 
It is noted that accurate frequency approximates are obtained with the cubic splines with 
a maximum relative error (over the reported results) of 0.3% occurring for the frequency 
2255.52 Hz. Due to the limited accuracy of the linear splines, significantly more basis functions 
must be used to obtain accurate approximations, and a 1.25% error at 2255.52 Hz remains, 
even with N — 32 functions. 

For n = 1, cross sections at 6 = 0 of corresponding modes approximated using M = 0, N — 
16 cubic splines are plotted in Figure 2. The coupling between the axial/radial components 
is readily verified for the 406.8 Hz and 603.83 Hz modes while the plot of the 266.05 Hz 
mode illustrates that for M = 0, the torsional mode is uncoupled (examples in Section 5.1.3 
illustrate that all three modes are coupled for M > 1). 

It should be noted that for each n , all three frequencies and modes are automatically 
yielded by the approximation method whereas (3.3) and (3.7) must be employed to obtain 
analytic expressions for the axial/radial and torsional modes, respectively. In this aspect, the 
approximation method facilitates the calculation of a complete modal set. 


n 


Analytic Frequencies 


Cubic Spline Approx 


Linear Spline Approx 


1 

2 

3 

4 

5 


266.05 


406.80 603.83 


531.53 
541.03 

543.53 
544.60 


532.11 


798.16 


1064.22 


1330.27 


266.06 


406.80 603.83 


266.19 


406.99 603.87 


| 924.29 

531.53 

532.11 

924.29 

531.71 

533.18 

925.95 

1362.11 

541.03 

■Em 

1362.19 

541.25 

801.79 

1368.07 

| 1807.86 

543.58 

1064.82 

1808.81 

543.85 

1072.81 

1822.21 

| 2255.52 

544.85 

1333.95 

2262.44 

545.06 

1347.10 

2283.78 


Table 3. Analytic frequencies and approximate values obtained using N = 8 cubic splines 
and N = 32 linear splines. Frequencies of the torsional modes are boxed. 
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Figure 2. Approximate axisymmetric modes U^{x, 0), V^(x, 0) and W^(x, 0) corresponding 
to the 544.60, 1330.27 and 2255.52 Hz natural frequencies. 
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5.1.2 m > l,n = 0 Purely Extensional Modes 

A second case in which analytic expressions for the shell frequencies and modes are easily ex- 
pressed occurs when purely extensional modes are present. The first three analytic frequencies 
given by (3.6) are compared in Table 4 with approximate values obtained with N = 8 cubic 
splines. Due to the accuracy of the method, the approximate frequencies agree with analytic 
values to within at least two decimal places in all three cases. Corresponding modal plots are 
given in Figure 3. The numerical plots in the left column illustrate the two approximate solu- 
tions, for each frequency, obtained by solving the matrix eigenvalue problem (4.13) (again, the 
full solution set is automatically obtained by the approximation method). The corresponding 
analytic modes given by (3.3) and (3.7) are depicted in the right column of the figure. With 
n = 0, it is seen from (3.3) and (3.7) that the analytic expressions for the longitudinal modes 
are U(x,6) = Ai cos(m0) and U(x,0) = A 2 sin(m0) with zero displacement expected in the 
circumferential and transverse components. This longitudinal behavior is noted for both the 
analytic and approximate modes in Figure 3. 



Analtyic (Hz) Galerkin (Hz) 

m = 1, n = 0 
m = 2, n — 0 
m — 3, n = 0 

338.75 338.75 

677.50 677.50 

1016.25 1016.25 


Table 4. Analytic and approximate frequencies for purely extensional modes obtained using 
N — 8 cubic splines. 

5.1.3 m > l,n > 1 General Shell Modes 

For the general case with m > l,n > 1, the axial, radial and torsional modes are all coupled 
with analytic values for the natural frequencies determined by (3.6) where 0^ n; is one of the 
three solutions to the cubic equation (3.5). These analytic values are compared in Table 5 
with approximate frequencies obtained using the cubic spline basis with N = 8, M = 3. For 
comparison sake, frequencies of the axisymmetric (m = 0) and purely extensional (n = 0) are 
also included in this table. The accuracy of the cubic spline discretization leads to approximate 
frequencies having relative errors less than 0.5% for the reported values (the largest relative 
error for the reported values occurs in the second frequency obtained with m = 3, n = 5). The 
convergence of the method and accuracy obtained using the cubic splines is further illustrated 
by the frequency results in Table 6 where approximates obtained with N = 8 and N = 16 are 
compared. A check of the relative errors shows that the method is converging more quickly 
than the expected G(h*) rate. 

The modes for this general case are again obtained from (3.3) and (3.7). The coupling 
between the axial, radial and torsional modes can be seen in Figure 4 where the 873.61 Hz 
mode is plotted along the axial line l = {(x,6)\6 = 7r/4, 0 < x < 12}. 

We reiterate that the full set of frequencies and modes can be approximated using the 
Fourier-Galerkin method by simply solving the matrix eigenvalue problem (4.13). Moreover, 
this technique can be directly extended to the problem of calculating frequencies and modes for 
shells with other boundary conditions by suitably modifying the spline basis. This is illustrated 
in the next section where the modal analysis of a shell with clamped ends is considered. 
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338.75 Hz 


677.50 Hz 


1016.25 Hz 


Figure 3. Purely extensional modes in the longitudinal direction. For N = 8 cubic splines 
numerical values of along the curve £ = {(x, 0)|O < 9 < 27 r, x = 0} are plotted in left column 
Corresponding shell modes given by (3.3) and (3.7) with n = 0 are depicted in right column 














Analytic Frequencies 

Galerkin Approximates 

m = 0 

7ft — 0 

ft = 1 

266.05 

406.80 

603.83 

ft = 1 

266.06 

406.80 

603.83 

CM 

II 

e 

531.53 

532.11 

924.29 

ft = 2 

531.53 

532.11 

924.29 

n ~ 3 

541.03 

798.16 

1362.11 

n — S 

541.03 

798.22 

1362.19 

n — 4 

543.53 

1064.22 

1807.86 

ft — 4 

543.58 

1064.82 

1808.81 

n = 5 

544.60 

1330.27 

2255.52 

ft = 5 

544.85 

1333.95 

2262.44 

m = 1 

7ft — 1 

n = 0 


338.75 


lEsa 


338.75 


n— 1 

147.09 

508.59 

873.61 

ft = 1 

147.09 

508.59 

873.61 

ft = 2 

327.96 

714.55 

1115.55 

n — 2 

327.96 

714.44 

1115.55 

ft — 3 

434.81 

909.78 

1486.93 

ft = 3 

434.83 

909.85 

1487.00 

n = 4 

483.84 

1137.35 

1900.34 

ft = 4 

483.92 

1138.09 

1901.13 

3 

II 

Cn 

507.08 

1383.62 

2329.09 

ft = 5 

507.48 

1388.30 

2335.07 


m 

= 2 



7ft 

= 2 


O 

II 

C 


677.50 


ft = 0 


677.50 


ft = 1 

64.56 

755.92 

1340.06 

ft = 1 

64.56 

755.92 

1340.06 

CM 

II 

e 

187.34 

917.53 

1521.01 

n = 2 

187.35 

917.53 

1521.01 

n = 3 

296.26 

1100.21 

1804.82 

n = 3 

296.28 

1100.29 

1804.85 

ft = 4 

373.46 

1300.70 

2153.02 

ft = 4 

373.61 

1301.67 

2153.52 

ft = 5 

423.98 

1519.53 

2536.80 

iO 

II 

c 

424.54 

1521.01 

2540.73 


m 

= 3 



7ft 

= 3 


ft = 0 


1016.25 


ft = 0 


1016.25 


ft = 1 

33.48 

1061.88 

1858.96 

ft = 1 

33.48 

1061.88 

1858.96 

n = 2 

111.09 

1177.72 

2001.04 

n = 2 

111.09 

1177.73 

2001.04 

n = 3 

198.51 

1332.12 

2225.92 

n — 3 

198.54 

1332.21 

2225.94 

ft = 4 

275.56 

1509.99 

2514.74 

ft = 4 

275.76 

1511.08 

2515.01 

ft = 5 

336.44 

1706.31 

2848.26 

ft = 5 

337.29 

1714.42 

2850.19 


Table 5. Frequencies approximated with M = 3 and N = 8 cubic basis functions. 


22 









m = 0, 

71—5 


N = 8 

544.85 

1333.95 

2262.44 

A = 16 

544.60 

1330.30 

2255.57 

Analytic 

544.60 

1330.27 

2255.52 


m = 1, 

71 = 5 


A = 8 

507.48 

1388.30 

2335.07 

N = 16 

507.08 

1383.66 

2329.14 

Analytic 

507.08 

1383.62 

2329.09 


m = 2, 

71—5 


N = 8 

424.54 

1521.01 

2540.73 

A = 16 

423.99 

1519.56 

2536.84 

Analytic 

423.98 

1519.53 

2536.80 


m = 3, 

77 — 5 


A = 8 

337.27 

1714.42 

2850.19 

A = 16 

336.45 

1706.35 

2848.28 

Analytic 

336.44 

1706.31 

2848.26 


Table 6. Frequencies approximated with M = 3 and N = 8,16 cubic basis functions. 



Figure 4. Approximate 873.61 Hz modes along the line £ = {(x,9)\0 = 7t/4,0 < x < 12}. 
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5.2 Modal Example: Fixed-Edge Shell 

The lack of analytic expressions for the frequencies of a clamped-end shell necessitates their 
numerical approximation. Like the simply supported case, the frequencies are approximated 
by solving the matrix eigenvalue problem (4.13) with basis functions constructed so as to 
satisfy the fixed-edge boundary conditions (2.6). Axisymmetric frequencies obtained using a 
Fourier limit of M = 0 and x-axis basis number of N = 16 are summarized in Table 7. A com- 
parison with the axisymmetric frequencies of the simply supported shell (see Table 3) shows 
that only the torsional frequencies match (the torsional frequencies are boxed to facilitate 
comparison). 


Approximate Frequencies 

266.05 

410.10 

532.11 

533.22 

541.76 

543.99 

544.96 

545.56 

546.07 

546.69 

547.63 

549.25 

552.18 

557.39 

565.57 

571.17 

575.20 

599.46 

798.16 

921.82 

1064.22 

1330.30 

1360.78 

1596.45 

1806.93 

1862.80 

2129.72 

2254.85 

2397.95 

2669.00 

2703.58 

2945.38 

3153.02 

3230.49 

3527.16 

3603.61 

3833.08 

4056.54 

4119.08 

4514.29 

4981.08 

5167.03 

5172.22 

5462.62 

5963.62 

6480.17 

8742.91 

6962.94 

8734.17 


Table 7 . Cubic spline approximates of axisymmetric frequencies obtained with M = 0 and 
N = 16 for the clamped-end shell. 

While the torsional modes for the fixed-edge shell match those for a simply-supported 
shell when m = 0 (see [16]), the general modes for the two shells differ due to boundary 
conditions. To illustrate, fixed-edge and simply-supported shell modes having similar fre- 
quencies and wavenumbers are compared in Figure 5. Specifically, the 545.56 Hz mode for 
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the clamped-edge shell is compared with the corresponding 545.66 Hz mode for the simply- 
supported shell. The plots show that while the qualitative mode shapes are somewhat similar 
for the two boundary conditions, the clamped-edge shell does not exhibit the trigonometric 
modal relations which characterize the simply-supported shell. This illustrates the neces- 
sity of obtaining approximations consistent with the physical application if employing modal 
expansions for simulations, damage detection or control. 



Figure 5. Axisymmetric modes U^(x, 0), V^(x, 0) and W^(x, 0) corresponding to the fixed- 
edge 545.56 Hz frequency and simply supported 545.66 Hz frequency. 

5.3 Forced Shell Examples 

A second means of testing the capabilities of the approximation method is through the consid- 
eration of steady state and time-dependent shell responses to given force and moment inputs. 
In the first examples of this section, known true solutions are used to calculate corresponding 
force input. By comparing the resulting approximate solutions with the original true solu- 
tions, the expected 0(h 2 x ) and 0(h x ) convergence rates for the linear and cubic splines can 
be verified. In the final example, a triangular input to the patch terms (modeling a voltage to 
spike to the patches) provides a broadband model response. The resulting frequencies are then 
compared with those obtained by solving the matrix eigenvalue problem (4.13) to illustrate 
the consistency among techniques. 
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5.3.1 Steady State Solution — Linear and Cubic Spline Approximations 

In this example, steady state, 0-dependent solutions satisfying the clamped edge boundary 
condition are considered. The true displacements were taken to be 

u(x,0) = sin(27rx/^) sin(0) 

v(x,0) = sin(47ra;/^)sin(20) ( 5 . 1 ) 

w(x,0) = [cos(27r;r/T) — l]cos(30) , 

and hence they satisfy the fixed-edge conditions (2.6). Plots of those displacements are given in 
Figure 6 to illustrate that significant bending deformations are present. The forcing functions 
were obtained by substituting the true solutions into the strong form (2.1) of the equations 
and solving for the functions q x ,q$ and q n . 

Since the forces are time-independent, the coefficients = [ui, • • • , ujq u , , • • • , u A ' . 
wi for the approximate displacements u N (0,x),v N (0,x) and w N (0*x) of (4.1) are 
obtained by solving the matrix system 


K% 0* = F* , 

where and F^ r are defined in (4.9) and (4.10), respectively. The specific construction 
and dimensions of Kg and F'^ depend upon the bases used in forming u N ,v N and w N (see 
Table 1). 

When approximating the solutions to (5.1), both the linear and cubic axial basis functions 
summarized in Table 1 were considered. For both cases, the maximum absolute errors over 
a uniform 25 x 25 grid in (x, 0) are summarized in Table 8 while ratios of the errors are 
given in Table 9 (the absolute errors provide the same information as relative errors since the 
magnitudes of the solutions are on the order of one). For both sets of approximations, equal 
axial index limits were assumed and are denoted by N u = N v = N w = N. Due to the nature 

of the forcing function, the Fourier limit M — 3 was sufficient for resolving the circumferential 
behavior. 

The error results in Table 8 illustrate the rapid convergence of the method with the fully 
cubic spline basis while the ratios in Table 9 demonstrate the respective 0(h *) and 0(h 2 ) 
convergence rates of the cubic and linear bases, respectively (the multiplying constant causes 
the ratios in the linear case to be slightly less than the expected value of four). It should be 
noted that due to the coupling between the longitudinal, circumferential and transverse dis- 
placements in the modeling equations, the accuracy of the transverse approximate is reduced 
to 0{h 2 x ) when linear splines were used to approximate u and v, even though cubic splines 
were used to approximate w. 

The approximate longitudinal displacements obtained with N = 4,8,16 and 32 linear 
axial basis functions are plotted in Figure 7. Comparison with the true solution in Figure 6 
illustrates that with N = 16, the approximation method has captured the qualitative behavior 
of the solution. Corresponding plots of the circumferential and transverse displacements, 
plotted in Figures 14 and 15, can be compared with the true circumferential and transverse 
solutions in Figure 6. 
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Fully Cubic Spline Axial Basis 


N = 4 


N = 8 


N = 16 


Mixed Linear/Cubic Axial Basis 


N = 4 


N = 8 


N = 16 


N = 32 


u 

V 

w 


1.663 - 1 
5.687 - 1 
3.136 - 1 


2.003 - 3 
1.501 - 2 
4.075 - 3 


3.936 - 5 
6.235 - 4 
5.552 - 5 


5.010-1 
1.580 + 0 
1.864 + 0 


1.266 -1 
2.962 - 1 
4.903 - 1 


3.277 - 2 
8.313 - 2 
1.236 - 1 


7.751 -3 
2.197-2 
3.784 - 2 


Table 8. Absolute errors in the longitudinal u, circumferential v , and transverse displace- 
ments w with cubic and linear spline axial basis functions. 


Cubic Basis 


N= 4 N=8 

K=8 N = 16 


Linear/Cubic Basis 

jV-4 N=S N=16 

N = 8 7V=16 N = 32 


V 


of 


w 


83.052 50.888 
37.899 24.069 
76.961 73.395 


3.956 3.863 4.228 
5.336 3.563 3.785 
3.801 3.966 3.268 




















Appro x*n»t* g, N-4 


Appfotfmow u. N-6 



Wmm 


mkmm 


Approxknolo u, N-t« 


ApproxJm*» u, N-32 


Figure 7. Approximate longitudinal displacement u with N = 4,8, 16,32 


Apprortmtto v, N-* 


Approtfrruto v'.f'A.tS 


ApproxJmat# v, N-32 


Thou 0 0 


Thou 0 0 


Figure 8. Approximate circumferential displacement v with N = 4, 8, 16 32 





AppfOJdrrwi* w, N-4 


Approximate w, N-8 




Approximate w, N»1 6 



Approximate w, N-32 



Figure 9. Approximate transverse displacement w with N = 4, 8, 16, 32. 


5.3.2 Time-Dependent Solution — Cubic Spline Approximation 

To illustrate the approximation of the dynamic shell model (including mass, stiffness, damping 
and force contributions), time- dependent forcing functions and solutions are considered in the 
final two examples. Cubic spline axial basis functions are used for discretizing the longitudinal, 
circumferential and tangential displacements in both cases. 

For this example, the true solutions were taken to be 

u(x , 9) = t 2 sin(27 tx/£) sin(0) 
v(x, 9) = i 1 sin(47 tx/£) sin(20) 
w(x,9) = t 2 [cos(27tx/£) — l]cos(30) 

(compare with (5.1)). Solving for q x ,<}9 and q n in the strong form of the equations (2.1) 
yields the forcing functions. Since the forcing terms are time- dependent, discretization of the 
modeling equations yields the matrix equation (4.12) in R n ,N = 2 Af. As detailed in [16], 
the resulting ODE system is moderately stiff and a variable-stepsize, variable-order stiff ODE 
routine was used to numerically integrate the system. 

The absolute errors in the longitudinal, circumferential and transverse displacements at 
time T — 1.5 sec are reported in Table 10 while error ratios as N was doubled are reported in 
Table 11. A comparison between the errors in Table 10 and the steady state errors in Table 8 
indicates that they are of the same magnitude while comparison between Tables 11 and 9 
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indicates nearly identical error ratios. Both Table 10 and 11 comfirm the 0(h *) convergence 
rate of the method in the damped, time-dependent problem. Hence with the fully cubic spline 
axial basis, the accuracy contributes to the efficiency of the method through reduced system 
sizes. This significantly speeds the time marching which is advantageous in simulations, 
parameter estimation and control applications. 



II 

N = 8 

N = 16 

u 

3.742 - 01 

4.506 - 03 

8.855 - 05 

V 

1.280 + 00 

3.376 - 02 

1.403 - 03 

w 

7.056 - 01 

9.169-03 

1.249 - 04 


Table 10. Absolute errors of the longitudinal u, circumferential v, and transverse displace- 
ments w at time T = 1.5 seconds. 



N = 4 

N = 8 

N = 8 

N = 16 

u 

83.052 

50.888 

V 

37.899 

24.069 

w 

76.961 

73.400 


Table 11. Ratios of the absolute errors at time T = 1.5 seconds. 


5.3.3 Time-Dependent Solution — Patch Input 

In the final example, the external inputs to the patch were assumed to come only from moments 
and forces generated by a single patch pair located between x x = 4.5 in, x 2 = 7.5 in and 
Q\ = O,0 2 = tt/ 6. Material and patch parameters are summarized in Table 12. The reported 
patch parameters are consistent with those of patches currently being used in the Acoustics 
Division, NASA Langley Research Center. 

For this case, the force vector F N (t ) of (4.10) has the components 


n x 2 < 97 ? 

R(N x ) pe (t)^ dx dd 
l ox 

(■N 0 ) pe (t )-gg - dx dd 

f $2 r r 

wh—i L 




dx d 6 


where ( N x ) pe , ( Ne ) pe , (M r ) pe and ( M$) pe are defined in (2.4) with i = 1. It is noted that due to 
the characteristic functions, the integrals reduce to only that region covered by the patches. 
This is in accordance with the physical fact that moments and forces are generated only in 
patch regions. 
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The voltage Vj(t) to the outer patch (see (2.5)) was taken to have a triangular shape 
with a magnitude of 95 V and a temporal duration of 1 x 10 -4 seconds. The voltage Vi(t) to 
the inner patch was taken to be V^t) = This input models an out-of-phase voltage 

spike to the patches and serves to generate a broadband shell response. The input of such 
voltages to patches is a current technique for eliciting shell responses when performing system 
identification or damage detection. 

The clamped-edge boundary condition (2.6) was assumed with compatible cubic basis 
functions in the numerical approximation. The index limits N = 8 and M = 3 were sufficient 
for resolving the shell dynamics in the previous examples, and those limits were also used to 
obtain the results here. The frequencies and displacements for the point (x,0) = (7.83,0.89) 
are plotted in Figures 10-12 (this point lies just outside the (x, 0) = (7.5,7 r/6) corner of the 
patch). Since out-of-phase simulated voltages were input to the inner and outer patches, the 
input was primarily in the form of external moments, with only higher order terms contributing 
to the in-plane forces (see the expressions of (2.5)). The small magnitudes of the displacements 
results from the use of a single actuating patch pair with a short duration voltage spike. Note 
that all computations were performed with double precision accuracy. 

In Table 13, the natural frequencies obtained through this dynamic patch excitation are 
compared with those obtained for the undamped shell by solving the matrix eigenvalue prob- 
lem (4.13). As expected, the two sets of frequencies are nearly identical since the same dis- 
cretization limits were used in both cases. The slight differences are due to the Kelvin-Voigt 
damping and sample rates in the dynamic simulation. 


h — 0.01 in 
p = 0.283 lb/in 3 
v = 0.3 
£ = 12.0 in 

R = 3.00 in 
E = 3x 10 7 lb/in 2 
cd — 15.099 lb • in • s 

Xi = 4.5 in 
0 ! =0 

E pe = 9.137xl0 6 lb/in 2 
h pe = 0.001 in 

X 2 — 7.5 in 
0 2 = 7 r/6 
I/ pe = 0.31 

d 3 1 = 7.4 xlO -9 in/V 


Table 12. Material and piezoceramic patch constants. 
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Figure 11. Circumferential displacements and frequencies for the damped shell. 
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Figure 12. Transverse displacements and frequencies for the damped shell. 


6 Concluding Remarks 

In this paper, a spline-based Galerkin method appropriate for discretizing cylindrical shell 
models has been presented. The emphasis throughout the discussion was on the development 
of a method which is sufficiently flexible for a variety of models including those which arise 
in smart material applications. Galerkin expansions, in terms of basis functions composed 
of modified splines in the axial direction and Fourier components circumferentially, provided 
this flexibility. Discretization in this manner also provided an accurate and efficient means 
of approximating static and dynamic shell deformations as well as natural frequencies and 
modes for a variety of boundary conditions. 

As demonstrated through numerical examples, the method yielded 0(h x ) convergence for 
physically realistic values of the shell thickness h when linear splines were used as axial basis 
elements for u,u, and 0(h x ) accuracy when all three components u,v,w were approximated 
by cubic spline expansions. The fourth order accuracy attained with cubic splines contributes 
to the method’s efficiency by reducing the number of basis elements and the hence system 
sizes necessary for resolving the shell behavior (in many of the examples, N = 8 cubic spline 
basis functions were sufficient for resolving the displacements u,v,w). Unlike the splines, the 
Fourier expansions in 6 form a truncated approximate in which all Fourier components with 
index less than the limit index M are resolved. In applications, the number of excited Fourier 
modes is often quite small, thus permitting small limits M. 

Another advantage of the cubic spline basis is the fact that differentiability and smooth- 
ness properties are constructed in the basis as compared with many finite element methods in 
which these conditions are obtained through compatibility conditions at the nodes. This sim- 
plifies construction of system matrices since components are constructed simply through the 
quadrature of basis elements. This also reduces the susceptibility of the method to membrane 
and shear locking since there are no elements per se in which asymptotic incompatibilities 
can arise. Alternatively, the spline-based Galerkin method can be thought of as similar to a 
spectral method but with basis functions having finite support. 
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Initial tests have indicated no tendency toward membrane or shear locking when using 
the fully cubic spline axial bases with small partition size h x . This is consistent with finite 
element results reported in [25]. These locking considerations, in combination with the 0(h *) 
accuracy of the cubic splines, indicate that the fully cubic spline axial basis will be preferable 
to the mixed linear/ cubic basis in most applications. 

The method was also demonstrated to be flexible with regards to boundary conditions, 
and both clamped-end and simply-supported edge conditions were attained through minor 
modifications in the boundary splines. Free end conditions are natural boundary conditions 
in the weak model formulation and hence no basis modifications are necessary in that case. 
“Almost fixed” boundary conditions modeling slight edge rotations and inplane movement are 
also common in applications due to the difficulty in maintaining perfect boundary clamps. As 
discussed in [7], only minor modifications are necessary for adopting the splines for this case. 

Finally, the approximation method is directly applicable to shell models incorporating pas- 
sive (internal) and active (external) contributions from embedded or surface-mounted sensors 
and/or actuators. When derived from physical principles, the passive actuator/sensor con- 
tributions are manifested in the physical model parameters (e.g., density, stiffness, damping, 
Poisson ratio) and are incorporated in the resulting discrete system simply through piecewise 
integration over the shell domain. The external moments and forces are isolated to regions 
covered by actuators through characteristic functions in the model. In the weak form, this 
yields integrals over the actuator regions. While partitions should be aligned with actuator 
edges to obtain optimal accuracy, neither the basis nor inner products must be modified to ac- 
count for sensors and/or actuators. This is advantageous when approximating shell dynamics 
in smart material applications. 
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